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ABSTRACT 

We consider the nonaxisymmetric hnear theory of radially-stratified disks. We work 
in a shearing-sheet-like approximation, where the vertical structure of the disk is ne- 
glected, and develop equations for the evolution of a plane-wave perturbation comoving 
with the shear flow (a shearing wave, or "shwave"). We calculate a complete solu- 
tion set for compressive and incompressive short-wavelength perturbations in both the 
stratified and unstratified shearing-sheet models. We develop expressions for the late- 
time asymptotic evolution of an individual shwave as well as for the expectation value 
of the energy for an ensemble of shwaves that are initially distributed isotropically in 
A:-space. We find that: (i) incompressive, short-wavelength perturbations in the unstrat- 
ified shearing sheet exhibit transient growth and asymptotic decay, but the energy of an 
ensemble of such shwaves is constant with time (consistent with Afshordi, Mukhopad- 
hyay, &: Narayan 2004); (ii) short-wavelength compressive shwaves grow asymptotically 
in the unstratified shearing sheet, as does the energy of an ensemble of such shwaves; 
(iii) incompressive shwaves in the stratified shearing sheet have density and azimuthal 
velocity perturbations 5S, 5vy ~ t~^^ (for |Ri| ^ 1), where Ri = N^/^qQ)"^ is the 
Richardson number, is the square of the radial Brunt- Vaisala frequency and qQ is 
the effective shear rate; (iv) the energy of an ensemble of incompressive shwaves in the 
stratified shearing sheet behaves asymptotically as Rit^~^^^ for |Ri| <C 1. For Keplerian 
disks with modest radial gradients, |Ri| is expected to be <C 1, and there will therefore 
be weak growth in a single shwave for Ri < and near-linear growth in the energy of 
an ensemble of shwaves, independent of the sign of Ri. 

Subject headings: accretion, accretion disks, solar system: formation, galaxies: nuclei 



1. Introduction 

Angular momentum transport is central to the evolution of astrophysical disks. In many disks 
angular momentum is likely redistributed internally by magnetohydrodynamic (MHD) turbulence 
driven by the magnetorotational instability (MRI; see Balbus &: Hawley 1998). But in portions 
of disks around young, low-mass stars, in cataclysmic-variable disks in quiescence, and in X-ray 
transients in quiescence (Stone, Gammie, Balbus, & Hawley 2000; Gammie & Menou 1998; Menou 
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2000), disks may be composed of gas that is so neutral that the MRI fails. It is therefore of interest 
to understand if there are purely hydrodynamic mechanisms for driving turbulence and angular 
momentum transport in disks. 

The case for hydrodynamic angular momentum transport is not promising. Numerical experi- 
ments carried out under conditions similar to those under which the MRI produces ample angular 
momentum fluxes- local shearing-box models- show small or negative angular momentum fluxes 
when the magnetic field is turned off ( Hawley, Gammie, &; Balbus 1995; Hawley, Gammie, & Bal- 
bus 1996). Unstratified shearing-sheet models show decaying angular momentum flux and kinetic 
energy when nonlinearly perturbed, yet recover the well known, high Reynolds number nonlinear 
instability of plane Couette flow when the parameters of the model are set appropriately (Balbus, 
Hawley, & Stone 1996; see, however, the recent results by Umurhan & Regev 2004). Local models 
with unstable vertical stratification show overturning and the development of convective turbulence, 
but the mean angular momentum flux is small and of the wrong sign (Stone & Balbus 1996). 

Linear theory of global disk models has long indicated the presence of instabilities associated 
with reflecting boundaries or features in the flow (see e.g., Papaloizou &; Pringle 1984, 1985, 1987; 
Goldreich, Goodman, &; Narayan 1986; Goodman, Narayan, & Goldreich 1987; Narayan, Goldreich, 
&; Goodman 1987; Lovelace, Li, Colgate, & Nelson 1999; Li, Finn, Lovelace, Sz Colgate 2000). 
Numerical simulations of the nonlinear outcome of these instabilities suggest that they saturate at 
low levels and are turned off by modest accretion (Blaes 1987; Hawley 1991). One might guess that 
in the nonlinear outcome these instabilities will attempt to smooth out the features that give rise to 
them, much as convection tends to erase its parent entropy gradient. There are some suggestions, 
however, that such instabilities saturate into long-lived vortices, which may serve as obstructions 
in the flow that give rise to angular momentum transport (Li, Colgate, Wendroff, &; Liska 2001). 
We will consider this possibility in a later publication. 

Linear theory has yet to uncover a local instability of hydrodynamic disks that produces 
astrophysically-relevant angular momentum fluxes. Because of the absence of a complete set of 
modes in the shearing-sheet model, however, local linear stability is difficult to prove. Local non- 
linear stability may be impossible to prove. Comparison with laboratory Couette flow experiments 
is complicated by several factors, not least of which is the inevitable presence of solid radial bound- 
aries in the laboratory that have no analogue in astrophysical disks. 

Recently, however, Klahr &: Bodenheimer (2003) (hereafter KB03) have claimed to find a local 
hydrodynamic instability in global numerical simulations. The instability arises in a model with 
scale-free initial conditions (an equilibrium entropy profile that varies as a power-law in radius) 
and thus does not depend on sharp features in the flow. Klahr (2004) has performed a local linear 
stability analysis of a radially-stratified accretion disk in an effort to explain the numerical results 
obtained by KB03. The instability mechanism invoked is the phenomenon of transient amplifica- 
tion as a shearing wave goes from leading to trailing. This is the mechanism that operates for 
nonaxisymmetric shearing waves in a disk that is nearly unstable to the axisymmetric gravitational 
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instability (Goldreich & Lynden-Bell 1965; Julian k, Toomre 1966; Goldreich k, Tremaine 1978). 
It is the purpose of this work to clarify and extend the linear analysis of Klahr (2004). If this 
instability exists it could be important for the evolution of low-ionization disks. 

To isolate the cause for instabilities originally observed in global 3D simulations, KB03 perform 
both local and global 2D calculations in the (/>)-plane. The local simulations use a new set 
of boundary conditions termed the shearing-disk boundary conditions. The model is designed to 
simulate a local portion of the disk without neglecting global effects such as curvature and horizontal 
flow gradients. The boundary conditions, which are described in more detail in KB03, require the 
assumption of a power-law scaling for the mean values of each of the variables, as well as the 
assumption that the fluctuations in each variable are proportional to their mean values. The radial 
velocity component in the inner and outer four grid cells is damped by 5% each time step in order 
to remove artificial radial oscillations produced by the model. ^ 

The equilibrium profile for KB03's 2D runs was a constant surface density S with either 
a constant temperature T or a temperature profile T oc R"^. The constant-T runs showed no 
instability while those with varying T (and thus varying entropy) sustained turbulence and positive 
Reynolds stresses.^ The fiducial local simulations were run at a resolution of 64^, with a spatial 
domain of i? = 4 to 6 AU and = 30°. The unstable run was repeated at a resolution of 
128^, along with a run at twice the physical size of the fiducial runs. One global model (with 
non-refiecting outflow boundary conditions) was run at a resolution of 128^ with a spatial domain 
of = 1 to 10 AU and A(/> = 360°. All the runs yielded similar results, with the larger simulations 
producing vortices and power on large scales. 

KB03 have chosen the term "baroclinic instability" by way of analogy with the baroclinic in- 
stability that gives rise to weather patterns in the atmosphere of the Earth and other planets (see 
e.g. Pedlosky 1979).^ The analogy is somewhat misleading, however, since the baroclinic instability 
that arises in planetary contexts is due to a baroclinic equilibrium. In a planetary atmosphere, a 
baroclinically-unstable situation requires stratification in both the vertical and latitudinal direc- 
tions.^ The stratification in KB03 is only in the radial direction, and as a result the equilibrium 
is barotropic. It is the perturbations that are baroclinic; i.e., the disk is only baroclinic at linear 
order in the amplitude of a disturbance. 



^It is not surprising that shearing disk boundary conditions as implemented in KB03 produce features on the 
radial boundary, because the Coriolis parameter is discontinuous across the radial boundary. 

^Notice that with a constant E, the constant-T runs have no variation in any of the equilibrium variables, so it is 
not clear that the effects being observed in the 2D calculations are due to the presence of an entropy gradient rather 
than due simply to the presence of a pressure gradient. 

baroclinic flow is one in which surfaces of constant density are inclined with respect to surfaces of constant 
pressure. If these surfaces coincide, the flow is termed barotropic. 

^Contrary to the claim in Klahr (2004), the two-layer model (Pedlosky 1979) does not ignore the vertical structure; 
it simply considers the lowest-order vertical mode. 
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Cabot (1984) and Knobloch & Spruit (1986) have analyzed a thin disk with a barochnic 
equiUbrium state (with both vertical and radial gradients) . The latter find that due to the dominant 
effect of the Keplerian shear, the instability only occurs if the radial scale height is comparable to 
the vertical scale height, a condition which is unlikely to be astrophysically relevant. As pointed out 
in KB03, the salient feature that is common to their analysis and the classical baroclinic instability 
is an equilibrium entropy gradient in the horizontal direction. As we show in §2, however, an 
entropy gradient is not required in order for two-dimensional perturbations to be baroclinic; any 
horizontal stratification will do. 

The "global baroclinic instability" claimed by KB03 is thus analogous to the classical baroclinic 
instability in the sense that both have the potential to give rise to convection.^ When neglecting 
vertical structure, however, the situation in an accretion disk is more closely analogous to a shearing, 
stratified atmosphere, the stability of which is governed by the classical Richardson criterion (Miles 
1961; Chimonas 1970). The only additional physics in a disk is the presence of the Coriolis force. 
Most analyses of a shearing, stratified atmosphere, however, only consider stratification profiles that 
are stable to convection. The primary question that Klahr (2004) and this work are addressing, 
then, is whether or not the presence of shear stabilizes a stratified equilibrium that would be 
unstable in its absence. 

We begin in §2 by outlining the basic equations for a local model of a thin disk. §§3 and 4 
describe the local linear theory for nonaxisymmetric sinusoidal perturbations in unstratified and 
radially-stratified disks, respectively. We summarize and discuss the implications of our findings in 
S5. 



2. Basic Equations 

The effect of radial gradients on the local stability of a thin disk can be analyzed most simply 
in the two-dimensional shearing-sheet approximation. This is obtained by a rigorous expansion of 
the equations of motion in the ratio of the vertical scale height H to the local radius i?, followed by 
a vertical integration of the fluid equations. The basic equations that one obtains (e.g., Goldreich 
&; Tremaine 1978) are 

f + EV.. = 0, (1) 
dv V P 

— + + 2nxv- 2qn'^xx = 0, (2) 



dt S 

dlnS 



dt 



0, (3) 



^The classical baroclinic instability gives rise to a form of "sloping convection" (Houghton 2002) since the latitu- 
dinal entropy gradient is inclined with respect to the vertical buoyancy force. 
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where S and P are the two-dimensional density and pressure, S = PS"''' is monotonically related 
to the fluid entropy,^ v is the fluid velocity and d/dt is the Lagrangian derivative. The third 
and fourth terms in equation (2) represent the Coriolis and centrifugal forces in the local model 
expansion, where O is the local rotation frequency, x is the radial Cartesian coordinate and q is 
the shear parameter (equal to 1.5 for a disk with a Keplerian rotation profile). The gravitational 
potential of the central object is included as part of the centrifugal force term in the local-model 
expansion, and we ignore the self-gravity of the disk. 

It is worth emphasizing at this point that we have integrated out the vertical degrees of freedom 
in the model. We will later focus on perturbations with planar wavelengths that are small compared 
to a scale height, and these perturbations will be strongly influenced by the vertical structure of 
the disk. 

Equations (1) through (3) can be combined into a single equation governing the evolution of 
the potential vorticity: 

dfVxv + 2n\ d$ VSxVP 

' ^ - ■■ (4) 



dt\ S J dt S3 ■ 

In two dimensions, ^ has only one nonzero component and can therefore be regarded as a scalar. 
Equation (4) demonstrates that for P = P(S) (as in the case of a strictly adiabatic evolution with 
isentropic initial conditions), the potential vorticity of fluid elements is conserved. For P ^ P(S), 
however, the potential vorticity evolves with time. A barotropic equilibrium stratification can result 
in baroclinic perturbations that cause the potential vorticity to evolve at linear order. This can be 
seen by linearizing the scalar version of equation (4): 

dSe ^r.^j, z • (VSq X VdP - VPq X V(^S) 

— + vo • V(5C + Sv ■ V^o = — ^ ^3 -, (5) 

where we have dropped the term oc VSq x VPq. Notice that an entropy gradient is not required 
for the evolution of the perturbed potential vorticity. For Sq = PqSq = constant, equation (5) 
reduces to 

-^ + vo.VS^ + 6v.V^,= ^ (6) 

Potential vorticity is conserved only in the limit of zero stratification (Pq = constant) or adiabatic 
perturbations {6S = 0). 



3. Unstratified Shearing Sheet 

Our goal is to understand the effects of radial stratification, but we begin by developing the 
linear theory of the standard (unstratified) shearing sheet, in which the equilibrium density and 



®With the assumptions of vertical hydrostatic equihbrium and neghgible self-gravity, the effective two-dimensional 
adiabatic index can be shown to be 7 = (87313 — 1)/(73d + 1) (e.g. Goldreich, Goodman, & Narayan 1986). 
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pressure are assumed to be spatially constant. This will serve to establish notation and method of 
analysis and to highlight the changes introduced by radial stratification in the next section. 

Our analysis follows that of Goldreich &; Tremaine (1978) except for our neglect of self-gravity. 
The equilibrium consists of a uniform sheet with E = Eq = constant, P = Pq = constant, and 
vq = —qQxy. Wc consider nonaxisymmetric Eulerian perturbations about this equilibrium with 
space-time dependence 6{t)exjp{ikx{t)x + ikyy), where 

kx{t) = kxo + q^kyt (7) 

(with kxo and ky > constant) is required to allow for a spatial Fourier decomposition of the 
perturbation. We will refer to these perturbations as shearing waves, or with some trepidation, but 
more compactly, as shwaves. 



3.1. Linearized Equations 

To linear order in the perturbation amplitudes, the dynamical equations reduce to 

— + ikxSvx + iky5vy = 0, (8) 
Eo 

SP 

Svx — 2Qdvy + ikx— = 0, (9) 
So 



dvy + (2 - q)n6vx + iky— = 0, (10) 



SP 

5'P 



+ Cg{ikx6vx + ikySvy) = 0, (11) 
So 

where = 7P0/S0 is the square of the equilibrium sound speed and an over-dot denotes a time 
derivative. 

The above system of equations admits four lincarly-indcpcndcnt solutions. Two of these are 
the non- vortical shwaves (solutions for which the perturbed potential vorticity is zero), which in 
the absence of self-gravity can be solved for exactly. The remaining two solutions are the vortical 
shwaves. When ky ^ the latter reduce to the zero-frequency modes of the axisymmetric version of 
equations (8) through (11). One of these (the entropy mode) remains unchanged in nonaxisymmetry 
(in a frame comoving with the shear). There is thus only one nontrivial vortical shwave in the 
unstratified shearing sheet. 

In the limit of tightly- wound shwaves (|A;a;| S> ky), the nonvortical and vortical shwaves are 

compressive and incompressive, respectively. In the short- wavelength limit (Hky ^ 1, where 
H = Cs/ft is the vertical scale height), the compressive and incompressive solutions remain well 
separated at all times, but for Hky < 0(1) there is mixing between them near kx = as an 
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incompressive shwave shears from leading to trailing. ' With the understanding that the distinction 

between compressive shwaves and incompressive shwaves as separate solutions is not valid for all 
time when Hky < 0(1), we generally choose to employ these terms over the more general but less 
intuitive terms "non-vortical" and "vortical." 

Based upon the above considerations, it is convenient to study the vortical shwave in the 
short- wavelength, low-frequency {dt <C Cgky) limit. This is equivalent to working in the Boussinesq 
approximation,^ which in the unstratified shearing sheet amounts to assuming incompressible flow. 
In this limit, equation (8) is replaced with 

kxSvx + ky5Vy = 0. (12) 

This demonstrates the incompressive nature of the vortical shwave in the short- wavelength limit. 



3.2. Solutions 

In the unstratified shearing sheet, equation (5) for the perturbed potential vorticity can be 
integrated to give: 

ikxSvy — ikySvx ^ ST, 
oiu = '■ ?0T^ = constant, (13) 

^0 ^0 

where = (2 — g)0/So is the equilibrium potential vorticity and we have employed the subscript 
u to highlight the fact that the perturbed potential vorticity is only constant in the unstratified 
shearing sheet. To obtain the compressive-shwave solutions, we set the constant (5^„ to zero. 
Combining equations (10) and (13) with = 0, one obtains an expression for 5vxc in terms of 
5vyc and its derivative: 

e , _ C%ky5Vyc - ^oTpdVyc , . 

OVxc- t2y2 , „2p ' ^^^> 

where the subscript c indicates a compressive shwave. The associated density and pressure pertur- 
bations are 

_ SPc _ . ^oTpkxSVyc + kySVyc , s 

''^c - 72" - ^ t2y-2 , „2Z,2 

via equation (13). Reinserting equation (14) into equation (10), taking one time derivative and 
replacing 5P via equation (11), we obtain the following remarkably simple equation: 

6vyc + (c^fe2 + K^) SVyc = 0, (16) 

where k"^ = k^ + ky and = (2 — q)^^^ is the epicyclic frequency. Changing to the dimensionless 
dependent variable 

'Icgky I kxo \ . 2csky 



'^This was pointed out to us by J. Goodman. 

*We demonstrate this equivalence in the Appendix. 



-8- 



and defining 

2qncsky ' ^^^^ 

the equation governing 6vyc becomes 

(fSVyc 



+ [-T'-C) Svy, = 0. (19) 



This is the parabolic cylinder equation (e.g. Abramowitz &; Stegun 1972), the solutions of which 
are parabolic cylinder functions. One representation of the general solution is 



6vyc = e 2 



i rp2 



c^M{--'-C,-,'-tA +c2Tm(--'-C,-,'-T^ 
U 2 ' 2' 2 7 V4 2 ' 2' 2 



(20) 



where ci and C2 are constants of integration and M is a confluent hypergeometric function. This 
completely specifies the compressive solutions for the unstratified shearing sheet, for any value of 
ky. 

Equation (19) has been analyzed in detail by Narayan, Goldreich, Sz Goodman (1987); their 
modal analysis yields the analogue of equation (19) in radial-position space rather than in the 
radial-wavcmimbcr (kx = kyT) space that forms the natural basis for our shwave analysis. One 
way of seeing the correspondence between the modes and shwaves is to take the Fourier transform 
of the asymptotic form of the solution. Appropriate linear combinations of the solutions given in 
equation (20) have the following asymptotic time dependence for r S> 1: 

^""yc ^ y^^^P {^^^^ V^*^^^ (^^ / (^skxdtj . (21) 

The Fourier transform of the above expression, evaluated by the method of stationary phase for 
Hky ^ 1, yields 

SvyciX)cx^e^p(±'-X^y (22) 

which is equivalent to the expressions given for the modes analyzed by Narayan, Goldreich, & Good- 
man (1987), in which the dimensionless spatial variable (with zero frequency, so that corotation is 
at a; = 0) is defined as 

X ^ J^x. (23) 

y Cs 

To obtain the incompressive shwave, we use the condition of incompressibility (equation (12)) 
to write 5vy in terms of 5vx, and then combine the dynamical equations (9) and (10) to eliminate 
5P. The incompressive shwave is given by: 

5vxi = Svxio-^, (24) 



-9- 



Svyi = -^Sv^i, (25) 

Ky 

' +2{q- l)n^^ , (26) 



So 7P0 icgky \ky c 

dicates an incompressive sh 
6vxi ai t = Oy This solution is uniformly valid for all time to leading order in [Hky)~^ <C 1. 



where the subscript i indicates an incompressive shwave, /cq = /c^q + ^"y 5vxiQ is the value of 



3.3. Energetics of the Incompressive Shwaves 

We define the kinetic energy in a single incompressive shwave as 

Eki = l^Svli + Svli) = l^oSvli'^ = l^oSvlio-^, (27) 

which peaks at kx = 0. This is not the only possible definition for the energy associated with 
a shear-flow disturbance; see the discussion in Appendix A of Narayan, Goldreich, &; Goodman 
(1987). 

One can also define an amplification factor for an individual shwave, 

_ Ekijkx = 0) _ ,klo /„„x 

which indicates that an arbitrary amount of transient amplification in kinetic energy can be obtained 
as one increases the amount of swing for a leading shwave {kxo ^ —ky)- This is essentially the 

mechanism invoked by Chagelishvili, Zahn, Tevzadze, & Lominadze (2003), Umurhan & Regev 
(2004) and Afshordi, Mukhopadhyay, Sz Narayan (2004) to argue for the onset of turbulence in 
unmagnetized Keplerian disks. 

Because only a small subset of all Fourier components achieve large amplification (those with 
initial wavevector very nearly aligned with the radius vector), one must ask what amplification 
is achieved for an astrophysically relevant set of initial conditions containing a superposition of 
Fourier components. It is natural to draw such a set of Fourier components from a distribution 
that is isotropic, or nearly so, when ko is large. 

Consider, then, perturbing a disk with a random set of incompressive perturbations (initial 
velocities perpendicular to ko) drawn from an isotropic, Gaussian random field and asking how the 
expectation value for the kinetic energy associated with the perturbations evolves with time. The 
evolution of the expected energy density is given by the following integral: 

(E,) = L' I d'ko{E,,) =L'I d'kol^o{Svl,o)j^- (29) 



® Chagelishvili, Zahn, Tevzadze, & Lominadze (2003) obtained this solution by starting with the assumption of 
incompressibility. 
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where () indicates an average over an ensemble of initial conditions, the first equality follows from 
Parseval's theorem, the second equality follows from the incompressive shwave solution (24)- (26) 
and therefore applies only for koH 3> 1, and is a normalizing factor with units of length squared. 

For initial conditions that are isotropic in feg {^Vxio = Sij±{kQ,9) sinO, where {Sy'^^ko)) is the 
expectation value for the initial incompressive perturbation as a function of A;o and tan^ = ky/kxo), 
the integral becomes 

1 f /'^'^ 1 

{Ei) = -SoL^ / kodko{Svl{ko)) / d9 . . --. (30) 

2 J Jo sm"^ 6 + {qQtsmO + cos6)^ 

Changing integration variables to r = qi}t + cot 0, the angular integral becomes 

oo 2 

dT — = 27r, (31) 

-oo J- + T 

which is independent of time; hence 

{Ei) = {Ei{t = 0)) (32) 
and we do not expect the total energy in incompressive shwaves to evolve. 

Although this result may appear to depend in detail on the assumption of isotropy, one can 
show that it really only depends on {E^iit = 0)) being smooth near sin 6 = 0, i.e. that there should 
not be a concentration of power in nearly radial wavevectors. This can be seen from the following 
argument. If we relax the assumption of isotropy, the angular integral becomes 

''\e^J^<^ ^. (33) 

/o sin2^ + (gm sine + COS 0)2 ^ ^ 

For q^lt ^ 1 the above integrand is sharply peaked in the narrow regions around tan 9 = —1/ (qQt) <^ 
1 (i.e., sm0 ~ 0). One can perform a Taylor-series expansion of {6v'j_(kQ,6)) in these regions, and 
as long as {6v'j_(ko, 9)) itself is not sharply peaked it is well approximated as a constant. A modest 
relaxation of the assumption of isotropy, then, will result in an asymptotically constant value for 
the energy integral. 

Based upon this analysis, large amplification in an individual shwave does not in itself argue for 
a transition to turbulence due to transient growth. One must also demonstrate that a "natural" set 
of perturbations can extract energy from the background shear flow. In the case of the unstratified 
shearing sheet, the energy of a random set of incompressive perturbations remains constant with 
time. This is consistent with the results of Umurhan & Regev (2004), who see asymptotic decay 
in linear theory, because they work with a finite set of wavevectors, each of which must decay 
asymptotically. 



f 

Jo 



3.4. Energetics of the Compressive Shwaves 



Here we calculate the energy evolution of the compressive shwaves for comparison purposes. 
We will consider the evolution of short-wavelength compressive shwaves in which only the initial 
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velocity is perturbed, both for simplicity and for consistency with our calculation of the short- 
wavelength incompressive shwaves. As before, we will assume that the initial kinetic energy is 
distributed isotropically. 

Wc use the WKB solutions to equation (16) with Hky ^ l}^ With the initial density per- 
turbation set to zero (consistent with our assumption of only initial velocity perturbations), the 
uniformly- valid asymptotic solution to leading order in {Hky)~^ is given by 

Svyc = SvycoJ — cos{W - Wo), (34) 
Svxc = j^Svyc, (35) 

Ky 

= 3irHc> (36) 

where the WKB eikonal is given by 

W = J Cskdt = ^ j v/r+^dr = ^(ryr+^ + ln(r+yr+7^)), (37) 



with Wq being the value of at i = 0. 
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Using equation (35), the energy integral for the compressive shwaves in the short- wavelength 
limit is ^ 

{E,) =L^J d^ko{Ekc) =L^J d^ol^oiSvl)'^. (38) 
With initial velocities now parallel to ko (and again isotropic), this becomes 

{Ec) = ^J^oL^ J kodkoiSvjiko)) dO sin^ 9 + {q^t sin 9 + cose f cos^{W - Wq). (39) 
For qVLt S> 1, the angular integral is approximated by 

/ de I sin 01 (1 + cos(2W - 2Wq)) ~ 2qM + .1^1^ cos(cskoqnt^ - 7r/4), (40) 
Jo V Cgko 

where the second approximation comes from employing the method of stationary phase. ""^^ In the 
short- wavelength limit, then, 

{Ec{qnt » 1)) = 2qnt {Ec{t = 0)) . (41) 



^''These solutions are the short-wavelength, /lijft-frequency {dt ^ 0{csky)) limit of the full set of linear equations 
in the shearing sheet; see the Appendix. 

^^This is not the same WKB solution that is calculated in the tight- winding approximation by Goldreich & Tremaine 
(1978); in that case Csky/ k <C 1, the opposite limit to that which we are considering here. The two WKB solutions 
match for T 1 in the absence of self-gravity. We have verified the accuracy of this solution by comparing it to the 
exact solution with acceptable results, and it is valid to leading order for all time. 

^^The first approximation breaks down near sin^ = 0, but the contribution of these regions to the integral is 
negligible for qOi » 1, in contrast to the situation for incompressive shwaves. 
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Thus the kinetic energy of an initially isotropic distribution of compressive shwaves grows, presum- 
ably at the expense of the background shear flow. 

The fate of a single compressive shwave is to steepen into a weak shock train and then decay. 
The fate of the field of weak shocks generated by an ensemble of compressive shwaves is less clear, 
but the mere presence of weak shocks does not indicate a transition to turbulence. 

4. Radially-Stratified Shearing Sheet 

We now generalize our analysis to include the possibility that the background density and 
pressure varies with x; this stratification is required for the manifestation of a convective instability. 
In order to use the shwave formalism we must assume that the background varies on a scale 
L H <^ Rso that the local model expansion (e.g., the neglect of curvature terms in the equations 
of motion) is still valid. 

With this assumption the equilibrium condition becomes 

where a prime denotes an x-derivative. One can regard the background flow as providing an effective 
shear rate 

q^l = -v'q (43) 

that varies with x, in which case = — J'' q{s)ds Qy. 

Localized on this background flow we will consider a shearing wave with kyL ^ 1. That is, we 
will consider nonaxisymmetric short- wavelength Eulerian perturbations with spacetime dependence 
S{t) cxp(i kx{t, s)ds + ikyU + ikzz), where ky and kz are constants and 

kxit, x) = kxo + q{x)9,kyt. (44) 

It may not be immediately obvious that this is a valid expansion since the shwaves sit on top 
of a radially-varying background (see Toomre 1969 for a discussion of waves in a slowly-varying 
background). But this is an ordinary WKB expansion in disguise. To see this, one need only 
transform to "comoving" coordinates x' = x, y' = y + q{s)dsQt, t' = t (this procedure may be 
more familiar in a cosmological context; as Balbus (1988) has pointed out, this is possible for any 
flow in which the velocities depend linearly on the spatial coordinates). In this frame the time- 
dependent wavevector given above is transformed to a time-independent wavevector. The price paid 
for this is that dx d^' + qfttdy' , so new explicit time dependences appear on the right hand side 
of the perturbed equations of motion, and the perturbed variables no longer have time dependence 
exp(iajt'). Instead, we must solve an ODE for 5{t'). The y' dependence can be decomposed as 
e^p{ikyy'). The x' dependence can be treated via WKB, since the perturbation may be assumed 



(42) 
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to have the form W{ex', et') exjp{ik' ■ x'). This "nearly diagonahzes" the operator dx'- Thus we are 
considering the evolution of a wavepacket in comoving coordinates — a "shwavepacket" . 

For this procedure to be valid two conditions must be met. First the usual WKB condition 

must apply, kyL ^ 1. Second, the parameters of the flow that are "seen" by the shwavepacket must 
change little on the characteristic timescale for variation of which is for the incompressive 
shwaves. For solid body rotation (g = 0) the group velocity (derivable from equation (59), below) 
is \vg\ < Nx/k (for positive squared Brunt- Vaisala frequency N^, defined below; for < 
the waves grow in place), so the timescale for change of wave packet parameters in this case is 
L/\vg\ > kL/Nx S> It seems reasonable to anticipate similarly long timescales when shear is 

present. As a final check, we have verified directly, using a code based on the ZEUS code of Stone 
& Norman (1992), that a vortical shwavepacket in the stratified shearing sheet remains localized 
as it swings from leading to trailing. 



4.1. Linearized Equations 

To linear order in the perturbation amplitudes, the dynamical equations reduce to 

(5S Sv 

— + -y^ + ikxSvx + ikySvy + ik^Svz = 0, (45) 

^0 -t^E 



5P_ _ 



5vx - 2n5vy + ikx— - y^— = 0, (46) 



6P 

5vy + (2 - q)Q.5vx + iky— = 0, (47) 

Eo 

6P 

Sv, + ik,— = 0, (48) 

— -c2— + c2-^ = 0, (49) 
Eo Eo Ls 

where 

— = -^ = — + — = ^ + ^ (5Q) 
Lp 7P0 Lj:^ Ls So 7-50 ^ ^ 

define the equilibrium pressure, density and entropy scale heights. We have included the vertical 

component of the velocity in order to make contact with an axisymmetric convective instability 

that is present in two dimensions, after which we will set kz to zero. 

We will be mainly interested in the incompressive shwaves because the short-wavelength com- 
pressive shwaves are unchanged at leading order by stratification. We will therefore work solely in 
the Boussinesq approximation.^^ In addition to the assumption of incompressibility, this approxi- 
mation considers SP to be negligible in the entropy equation; pressure changes are determined by 



^We also drop the subscripts distinguishing between the compressive and incompressive shwaves. 
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whatever is required to maintain nearly incompressible flow. The original Boussinesq approxima- 
tion applies only to incompressible fluids. It was extended to compressible fluids by Jeffreys (1930) 
and Spiegel & Veronis (1960). We show in the Appendix that it is formally equivalent to taking the 
short-wavelength, low-frequency limit of the full set of linear equations. From this viewpoint, as- 
suming that HkySP/Po is of the same order as the other terms in the dynamical equations implies 
that SP/Pq ~ {Hky)~^S'S/'S,o, thus justifying its neglect in the entropy equation. We therefore 
replace equations (45) and (49) with 

kxSvx + kySvy + kz5vz = (51) 

and 

(5S 6vx 

— - = 0. (52) 
Using equations (47) and (52) and the time derivative of equation (51), one can express Svy 

5P .kx6vx + 2{q — l)QkySvx 



and 6P in terms of Svx and Svx- 



Yjq ky -\- k"^ 



(53) 



i-Qky + iq- 2)ki)nSvx - kxkySvx 
^""y = ■ (54) 

Eliminating SP in equation (46) via equation (53) gives 

P5vx + 2{q - l)nkxky5vx = {kl + k'^z){2n5vy + {cI/Lp)5T./Y.q), (55) 

where = k^ + k'^ + kl- Taking the time derivative of this equation and eliminating 5S and 5vy 
via equations (52) and (54), we obtain the following differential equation for 5vx- 

PSvx + Aqilkxky&vx + [kl {Nl + 2f^^) + k^ {N^ + k^)] 6vx = 0, (56) 

where = 2{2 — g)r2^ is the square of the effective epicyclic frequency and 

N'x = -j^ (57) 
is the square of the Brunt- Vaisala frequency in the radial direction. 



^''Notice that N^, q and /t^ are all functions of x and vary on a scale L H. 
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4.2. Comparison with Known Results 

Setting ky = in equation (56) yields the axisymmetric modes with the following dispersion 
relation (for d{t) oc e'^'^y. 

-' = TATrA^l + -')- (58) 



This is the origin of the H0iland stability criterion: the axisymmetric modes are stable for N'^ + k? > 
0. In the absence of rotation this reduces to the Schwarzschild stability criterion: A'^^ > is the 
necessary condition for stability. The effect of rotation is strongly stabilizing: if < —k^, as 
required for instability, then LsLp ~ H^; pressure and entropy must vary on radial scales of order 
the scale height for the disk to be H0iland unstable. 

Notice that effective epicyclic frequency only stabilizes modes with nonzero k^- The stability 
of nonaxisymmetric shwaves with kz = (as in the mid-plane of a thin disk) is the open question 
that this work is addressing. In this limit and in the absence of shear the Schwarzschild stability 
criterion is again recovered: with kz = and q = in equation (56) the dispersion relation becomes 

If there is a region of the disk where the effective shear is zero, a WKB normal-mode analysis will 
yield the above dispersion relation and there will be convective instability for ATj < 0. It appears 
from equation (56) that differential rotation provides a stabilizing influence for nonaxisymmetric 
shwaves just as rotation docs for the axisymmetric modes. Things are not as simple in nonax- 
isymmetry, however. The time dependence is no longer exponential, nor is it the same for all the 
perturbation variables. There is no clear cutoff between exponential and oscillatory behavior, so 
the question of flow stability becomes more subtle. 

As discussed in the introduction, the Boussinesq system of equations in the shearing-sheet 
model of a radially-stratified disk bear a close resemblance to the system of equations employed in 
analyses of a shearing, stratified atmosphere. A sufficient condition for stability in the latter case 
is that 

N'^ 1 

Ri = — % > - (60) 

everywhere in the flow, where Ri is the Richardson number, a measure of the relative importance 
of buoyancy and shear. This stability criterion was originally proved by Miles (1961) and Howard 
(1961) for incompressible fluids, and its extension to compressible fluids was demonstrated by 
Chimonas (1970). The stability criterion is based on a normal-mode analysis with rigid boundary 
conditions. Other than differences in notation (e.g., our radial coordinate corresponds to the vertical 
coordinate in a stratified atmosphere), the key differences in our system are: (i) the equilibrium 
pressure gradient in a disk is balanced by centrifugal forces rather than by gravity; (ii) the disk 
equations contain Coriolis force terms; (iii) most atmospheric analyses only consider an equilibrium 
that is convectively stable, whereas we are interested in an unstable stratification; (iv) we do not 
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employ boundary conditions in our analytic model since we are only interested in the possibility of 
a local instability. 

The lack of boundary conditions in our model makes the applicability of the standard Richard- 
son stability criterion in determining local stability somewhat dubious, since the lack of boundary 
conditions precludes the decomposition of linear disturbances into normal modes. The natural 
procedure for performing a local linear analysis in disks is to decompose the perturbations into 
shwaves, as we have done. 

Eliassen, H0iland, &: Riis (1953) consider both stable and unstable atmospheres and analyze an 
initial- value problem by decomposing the perturbations in time via Laplace transforms. For flow 
between two parallel walls, they find that an arbitrary initial disturbance behaves asymptotically 
as for -3/4 < Ri < 1/4, where 

a = Vl -4Ri, (61) 

which grows algebraically for Ri < 0. The disturbance grows exponentially only for Ri < —3/4. 
For a semi-infinite flow, the power-law behavior in time holds for —2 < Ri < 1/4, with exponential 
growth for Ri < —2. These results illustrate the importance of boundary conditions in determining 
stability. 

In the kz = limit that we are concerned with here, the correspondence between the disk and 
atmospheric models turns out to be exact in the shwave formalism. This is because the Coriolis 
force only appears in equation (56) via k^, which disappears when kz = 0. The equation describing 
the time evolution of shwaves in both a radially-stratified disk and a shearing, stratified atmosphere 
is thus 

Pdv^ + AqQk^kySv^ + {Nl + 2q^^^) 5v^ = 0. (62) 
We analyze the solutions to this equation in the following section. 



4.3. Solutions 

Changing time variables in equation (62) to f = kx/ky, the differential equation governing Sv^ 
becomes 

(1 + f2)4§^ + 4f^ + (Ri + 2)Svx = 0. (63) 

The solutions to equation (63) are hypergeometric functions. With the change of variables z = — f^, 
equation (63) becomes 

(P6vx l-5zdSvx Ri + 2 



^^Cited in Miles (1961). 
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The hypergeometric equation (Abramowitz &; Stegun 1972) 

z{l - z)'^ + [c - (a + 6 + l)z] ^ - ahSv, = (65) 

has as its two linearly independent solutions F{a, b; c; z) and z^~'^F(a — c + 1, 6 — c + 1; 2 — c; z). 
Comparison of equations (64) and (65) shows that a = (3 — a)/4, b= (3 + a)/4 and c = 1/2, where 
a is defined in equation (61). 

The general solution for 5vx is thus given by 

3 — a 3 + a 1 ~ rn f ^ ~ S + a 3 ,2 



Sv. = C,F -; -f^ j + f F (^^, -; -f , (66) 

where Ci and C2 are constants of integration representing the two degrees of freedom in our 
reduced system. These two degrees of freedom can be represented physically by the initial velocity 
and displacement of a perturbed fluid particle in the radial direction. The radial Lagrangian 
displacement is obtained from equation (66) by direct integration/^ 

^ /a ^2 p/l-a 1 + a 1 Ci + 

^^ = y^"^^* = -|j^m^(,^'^'2'-" J+^"^(,^'^'2'-" J- 

The solutions for the other perturbation variables can be obtained from equations (51), (52) and 
(53) with = 0: 

Svy = - fSvx, (68) 

(69) 

and 



SP _ jn 



(5S 






Ls 






^^Tf ( 


. Cs J 



dVx 
) — 



(70) 



It can be seen from the latter equation and the solution for Svx that SP/ Pq remains small compared 
to Svx/cg in the short- wavelength limit. This demonstrates the consistency of the Boussinesq 
approximation. 

The hypergeometric functions can be transformed to a form valid for large f (see Abramowitz 
&L Stegun 1972 equations 15.3.7 and 15.1.1). An equivalent form of the solution for |f| S> 1 is 

Svx = iCV, + sgn(f)C2F2) Irl^F 1 - ^; + 



(CiFs + sgn(T)C2F4) |t| 2 f{—^,—^;1 + -;-^) , (71) 



^®In our notation, a subscript x or y on the symbol ^ indicates a Lagrangian displacement, not a component of the 
potential vorticity, which is a scalar. 
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where sgn(f) is the arithmetic sign of f and the constants Vi are given by 

r(i)r(f) r(i)r(f) 
r(3±^)r(-i-^) ' - r (5±r.) r (i±^) ' 
y ^ rQ)r(-f) ^ r(i)r(-f) 

Expanding the above form of the solution for |f| 3> 1, we obtain 

Sv^ = {CiVi + sga{f)C2V2) \f\^ + (C1F3 + sgn(f)C2y4) \f\-^ + 0{f-^). (73) 



An equivalent form of $,x for |f| 3> 1 is 



C2X1 ,~^Cl^2^ ,.,H^„/3-a 1-a ^ a 1 



^+sgn(r)— — |r| 2 F [ ——, ——■,1 + ] , (74) 



gORi ^ ^ ' qn J ^ ^ V 4 ' 4 ' 2' f2 

where the constants Xi are given by 

- rQ)r(f) ^_ r(i)r(f) 

J- -n / 1 -Urv \ T-> / 1 -krv \ ' 



rmr(i±^) ' ^ - r (3±^) r (3+«) ' 
^ = r(^)r(-f) ^ r(|)r(-f) 

r(i-^)r(V) ' r(3^)r(3^)- ^''^ 

Expanding ^-^ for |f| S> 1 yields 



The dominant contribution for each perturbation variable at late times is thus 

5P oc Sva, t"^ , (77) 
S^cx^x^t^, (78) 

and 

a— 1 

6vy (X tSvx ~ . (79) 

This leads to one of our main conclusions: the density and y-vclocity perturbations will grow 
asymptotically for a > 1, i.e. Ri oc < 0.^^ For small Richardson number, however (as is 



^^Notice that this is the same time dependence obtained by Eliassen, H0iland, & Riis (1953) in a modal analysis; 
see the discussion surrounding equation (61). These power law time-dependences can be obtained more efficiently by 
solving the large-f limit of equation (66). 
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expected for a Keplerian disk with modest radial gradients), a ~ 1 — 2Ri and the asymptotic 
growth is extremely slow: 

~ 5vy - t-^\ (80) 

In the stratified shearing sheet, the right-hand side of equation (5) governing the evolution of 
the perturbed potential vorticity is no longer zero. The form of this equation for the incompressive 
shwaves is 

dt \ So I iLpTiQ 

The asymptotic time dependence of the perturbed potential vorticity can be obtained by integrating 
equation (81): 

S^^t"^^ t'-^^ (82) 

for f 3> 1 and |Ri| <C 1. As noted in §2, an entropy gradient is not required to generate vorticity. 
For iVj = 0, a = 1 and the perturbed potential vorticity grows linearly with time. The unstratified 
shearing sheet is recovered in the limit of zero stratification (1/Lp — 0), since in this limit equation 
(81) reduces to ^ = constant. 



4.4. Energetics of the Incompressive Shwaves 

For a physical interpretation of the incompressive shwaves in the stratified shearing sheet, we 
repeat the analysis of section 3.3 for the solution given in the previous section. For a complete 
description of the energy in this case, however, we must include the potential energy of a fluid 
element displaced in the radial direction. Following Miles (1961), an expression for the energy 
in the Boussinesq approximation is obtained by summing equation (46) multiplied by Sv^ and 
equation (47) multiplied by Svy. Replacing (5S/So by ^x/Ls via equation (69) results in the 
following expression for the energy evolution: 

^ " ^ (^^°'^"' + ^^oiVje^) = J^oSvJvy, (83) 

where 6v^ = dv^ + 5vy . The three terms in equation (83) can be identified as the kinetic energy, 
potential energy and Reynolds stress associated with an individual shwave. One may readily verify 
that the vortical shwaves (see equations (24)-(26)) in the unstratified shearing sheet {N"^ = 0) 
satisfy equation (83). 

The right hand side of equation (83) can be rewritten —f5v^ and individual trailing shwaves 
(f > 0) are therefore associated with a negative angular momentum flux. If the energy were positive 
definite this would require that individual shwaves always decay. But when A'^^ < (Ri < 0) the 
potential energy associated with a displacement is negative, so the energy Ek can be negative and 
a negative angular momentum flux is not enough to halt shwave growth. 
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Our next step is to write the constants of integration Ci and C2 in terms of the initial radial 
velocity and displacement of the shearing wave, 5vxo and ^xQ'- 



5vxi (fo) ixi (to) + Ri 5vx2 (to) ix2 (ro) 



-gORifaa;i(fo) ^xQ + Ri^x2(ro) 5vxo 

&Vx\{fQ) ^xAtq) + K\5Vx2{fQ) ^x2{fo) ' 



where fo = kxo/ky, Svxi is the hypergeometric function given by equation (66) with Ci = 1 and 
C2 = 0, and the other functions arc similarly defined. These expressions can be simplified by 
noticing that the denominator of Ci and C2 is the Wronskian of the differential equation for (x-^^ 



The Wronskian of this equation is 



yy = —p- t,xi t,x2 = exp 

ctr ar 



2r^ 



■dT 



1 + r2 J 1 + f2 ■ 
We further simplify the analysis by setting the initial displacement ^a;o to zero. 

With these simplifications, the solution given by equations (66) and (67) becomes 



8v^. -oxT^/l — a l + a 1 Tp(^~^ 3 + a 1 ^2 



5v 



'xG 



= (1 + ^0^) 



4 ' 4 '2' 



RifoF 



3 — a 3 + a 3 



4 ' 4 '2 



4 ' 4 '2' 
5 — a 5 + a 3 



4 ' 4 '2' 



(85) 



(86) 



(87) 



5v. 



'xO 



1 ^/3 — a3 + a3 oX^/l — al + al 9 



4 ' 4 '2' 



1 — a 1 + a 1 ^„/3 — a 3 + a 3 ^ 



(88) 



As in section 3.3, the energy integral for the incompressive perturbations is given by 

t2 I 



{Ei) = l^oL^J kodko{5vl{ko)) \esin^e (l + f^) (^^^ + iV, 



^ \ Sv. 



xO 



(89) 



for initial perturbations perpendicular to and isotropic in kg. Changing integration variables to 
f = qUt + cot 9, the angular integral becomes 



/OO |- 
df (1 + f2) {e^i(f - gOi)<^^xl(T) + Ri^a.2(f - qnt)dVx2{f)y + 
-00 



Ri {Cx2if - qnt)^xiif) - uif - m)^x2if)y 



(90) 



^®Based upon the relationship between a hypergeometric function and its derivatives, Svxi = d{^x2)/df and 
RiSvx2 = —d{ixi)/df. 
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where we have used the relation sin^ = (1 + Tq)"^. In the hmit of large qUt, the dominant 
contribution to the angular integral comes from the region < f < qO,t. This can be seen from the 
following argument. Using the expansions given by equations (73) and (76), we find the angular 
integrand is 

2|f(f-gJ^t)|"-i [{ViXi + sgn(f)sgn(f - qnt)mV2X2)^ + RiX2x|(sgn(f) - sgn(f - qnt))^] (91) 
for |f| S> 1 and |f — qClt\ S> 1. Using the relation T{n + 1) = nr{n), one can easily show that 

X2 = ^—Vi and V2 = ^—Xi. (92) 
a — 1 a + 1 

The integrand therefore simplifies to 

|f(f - qnt)\'^-%^Xf-^ [sgn(f) - sgn(f - qnt)f , (93) 
1 — a 

which is zero unless < f < qUt (for t > 0). For large qflt, therefore, the angular integral is 
approximately given by 

16Vf X2 ~oM«-l IQV^^X^ {fqntr^f ^ ^ ~ ^ q^t-u 



.6^2x2 fi^'t-'^ , ^ ^ ^ 16^.2x2 (fqnt)^ ^ / f \ 



(94) 



where 1 ^ ^ For S> i^, the above expression can be approximated by evaluating it at 
f = qilt, giving 

{E,{qm » 1)) l6V,'xf^^^±^^{qntf''-' ^ = 0)), (95) 

where we have used equation 15.3.7 in Abramowitz & Stegun (1972) to evaluate F(a, 6;c; 1).^^ 

Notice that there is no power- law growth in the perturbation energy for Ri > 1/4,^° consistent 
with the classical Richardson criterion (60) . In our analysis the energy decays with time for 2a — 
1 < 0, or Ri > 3/16. Thus the energy of an initial isotropic set of incompressive perturbations 
in a radially-stratified shearing sheet-model grows asymptotically (for Ri < 3/16), just like the 
compressive shwaves and unlike the incompressive shwaves in an unstratified shearing sheet, for 
which the energy is constant in time. 

The growth of an ensemble of incompressive shwaves in a stratified disk is not due to a Raylcigh- 
Taylor or convcctive type instability. There is asymptotic growth for < Ri < 3/16, and convective 
instability requires Ri < 0. One can also see this by examining the asymptotic energy for small 
values of |Ri|, such as would be expected for a Keplerian disk with modest radial gradients: 

{Eiiqnt > 1)) [2^2j^i + 0(Ri2)] {Ei{t = 0)). (96) 



^^We have numerically integrated the angular integral (90) and found this to be an excellent approximation at late 
times. 

^opor Ri > 1/4, a is imaginary and Re[t^"-^] = cos(2|a| Int). 
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Evidently for small values of Ri the near-linear growth in time of the energy is independent of the 
sign of Ri and therefore iVj.^^ 

5. Implications 

We have studied the nonaxisymmetric linear theory of a thin, radially-stratified disk. Our 
findings are: (i) incompressive, short-wavelength perturbations in the unstratified shearing sheet 
exhibit transient growth and asymptotic decay, but the energy of an ensemble of such shwaves is con- 
stant with time (consistent with Afshordi, Mukhopadhyay, & Narayan 2004) ; (ii) short- wavelength 
compressive shwaves grow asymptotically in the unstratified shearing sheet, as does the energy of 
an ensemble of such shwaves, which in the absence of any other dissipative effects (e.g., radiative 
damping) will result in a compressive shwave steepening into a train of weak shocks; (iii) incompres- 
sive shwaves in the stratified shearing sheet have density and azimuthal velocity perturbations 5T,, 
6vy ~ (for |Ri| < 1); (iv) incompressive shwaves in the stratified shearing sheet are associated 
with an angular momentum flux proportional to —kx/ky; leading shwaves therefore have positive 
angular momentum flux and trailing shwaves have negative angular momentum flux; (v) the energy 
of an ensemble of incompressive shwaves in the stratified shearing sheet behaves asymptotically as 
^i-4Ri £qj, ^ -p^j. Keplerian disks with modest radial gradients, |Ri| is expected to be ^ 1, 
and there will therefore be weak growth in a single shwave for Ri < and near-linear growth in the 
energy of an ensemble of shwaves, independent of the sign of Ri. 

Along the way we have found the following solutions: (i) an exact solution for non-vortical 
shwaves in the unstratified shearing sheet, equations (14), (15) and (20); (ii) a WKB-solution 
for the non-vortical, compressive shwaves in the short-wavelength, high-frequency limit, equations 
(34)- (36); (iii) a solution for incompressive shwaves in the unstratified shearing sheet valid in the 
short- wavelength, low-frequency limit, equations (24)-(26); (iv) a solution for incompressive shwaves 
in the radially-stratified shearing sheet (also valid in the short- wavelength, low- frequency limit), 
equations (66)- (70). 

Our results are summarized in Figure 1, which shows the regions of amplification and decay 
for shwaves in a stratified disk in the N^/^, q plane. 

The presence of power-law growth of incompressive shwaves in stratified disks opens the pos- 
sibility of a transition to turbulence as amplified shwaves enter the nonlinear regime. Any such 
transition would depend, however, on the nonlinear behavior of the disk after the shwaves break. 
It is far from clear that they would continue to grow. We will evaluate the nonlinear behavior of 
the disk in subsequent work. 

Our results are essentially in agreement with the numerical results presented by Klahr (2004) 



This asymptotic expression assumes Ri ^ 0. Notice that the energy at late times can have the opposite sign to 
the initial energy because the potential energy is negative for < 0. 
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(although the decay of high frequency, compressive shwaves seen in his Figures 8, 9, 10 is likely 
a numerical effect since compressive shwaves grow asymptotically in our analysis), that is, we 
find that arbitrarily large amplification factors can be obtained by starting with appropriate initial 
conditions. Our results, however, clarify the nature and asymptotic time dependence of the growth. 
Our results on the unstratified shearing sheet are also consistent with the results of Afshordi, 
Mukhopadhyay, k, Narayan (2004), who find that an isotropic ensemble of incompressive shwaves 
have fixed energy. 

This work was supported by NSF grant AST 00-03091 and a Drickamer Fellowship for BMJ. 
We thank Jeremy Goodman for his comments. 
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We demonstrate here that the Boussinesq approximation to the Unear perturbation equations 
is formally equivalent to a short- wavelength, low-frequency limit of the full set of linear equations. 
We perform the demonstration for the stratified shearing-sheet model since the standard shearing 
sheet is recovered in the limit of zero stratification. 

Combining equations (45) through (49) into a single equation for 6vx yields the following 
differential equation, fourth-order in time: 



where 



F4 



F^Svi^^ + ^3^4^) + ^2^4^) + FiSvi^^ + FoSv^ = 0, 

2 2(g+l)(g + 2) 



iK + k. 



1 - ^ 



kxLp 



F3 = -2qnkxky{k^ + ki)ll + j 



kxLp 



(1) 

(2) 
(3) 



F2 = clkl 



k;i\i + 



{k'y + ki) (fe^ + ki) 1 



kxLp 



kxLp 



+ k: 



, 2(g + l)(3g + 2) 



(4) 



Fi = 4:qnc%kl 



{k' + ki){l + 



kxLp, 



^klLl 

2(g + l)(g + 2) 



2qkxLp 
3k 



Fo = clkl [kliNl + 2cfn^) + kl{Nl + i^)] 



{ky + k^) 



klm 



kxLp 



2kim 

2i 



+ 



-I- k 



1 - 

SqkxLpJ _ 
2 2(g + l)(3^ + 2) 



k^H'^ 



(5) 
(6) 



The above expressions have been written to make the short-wavclcngth limit more apparent: 
all but the leading-order terms in brackets are proportional to factors of {k^Lp) '^ or {kxH)^^. 
Notice also that since one expects H/Lp <^ 1 for Keplerian disks with modest radial gradients, 



1 H 



kxLp kyHf Lp 



1 



kyHf^ 



(7) 



(for f 7^ 0) and therefore the short- wavelength limit is sufficient. One needs to be careful in taking 
this limit, however, f goes through zero as a shwave goes from leading to trailing. 

The approximation is rigorously valid only for f 7^ 0, but we have numerically integrated the full 
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set of linear equations (equations (45) through (49) with kz = 0) and found good agreement with 
the Boussinesq solutions described in §4.3 for all f at sufficiently short wavelengths.^^ 

With these assumptions in mind, to leading order in {Hky)~^ equation (1) becomes 

kjv^j'^ - 2qVlky5vf'^ + c%P5vP + 4qnciklky5v^^^ + 

c% [kl{N^ + 2fn^) + kliNl + R")] 5v^ = 0. (8) 

If we assume dt <^ Cgky, the two highest-order time derivatives are of lower order and can be 
neglected (thereby eliminating the compressive shwaves) and we have 

Pd'v:c + ^nk:ckySv^ [k^iN^ + 2q'^n'^) + A;2(Ar| + k^)] Sv^ = 0. (9) 

This is equivalent to equation (56). 

Notice also that the assumption dt ~ 0{csky) applied to equation (8) yields 

Svi^^ + clk^Svf^ = (10) 

to leading order in (Hky)^^. This equation is of the same form as the short-wavelength limit of 
equation (16) for the compressive shwaves in the unstratified shearing sheet, confirming our claim 
that short-wavelength compressive shwaves are unchanged at leading order by stratification. 



^^One must start with a set of initial conditions consistent with equations (46), (47), (51) and (52) in order to 
accurately track the incompressive-shwave solutions. In addition, suppression of the high-frequency compressive- 
shwave solutions near f = requires kyLp > 200, which for H/Lp = 0.1 implies Hky > 20. 
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Fig. 1. — A summary of analytic results for shwaves (shearing waves) in a stratified disk. The 
relevant parameters are the local shearing rate q = —^dlnQ.'^/dliir and the dimensionless Brunt- 
Vaisala frequency N'^/Vl?. The expected location of a thin, smooth disk is shown as a vertically 
extended ellipse near g = 0, iVj/O^ = 0. The far right region (shaded in the figure) is forbidden by 
the H0iland criterion. When g = shear is absent and a modal analysis is possible; instability is 
present for TVj < 0. Solitary shwaves with Ri = N^/{(ffP) < experience asymptotic power-law 
growth (oc t^^^ for small Ri); since each shwave grows the energy of an ensemble of shwaves docs 
as well. For < Ri < 3/16 solitary shwaves decay but the energy of an ensemble of shwaves grows 
as a power-law in time. For Ri > 3/16 both solitary shwaves and the energy of an ensemble of 
shwaves asymptotically decay. 



